#include<stdio.h>  
#include<stdlib.h>  
#include<math.h>  
#include<string.h>  
#include<time.h>
#define MAX_POINTS 400
#define EPOCHS 1000

typedef struct {  
     double x;  
     double y;  
} Point;  


void least_squares(const char* filename) 
{
    Point points[MAX_POINTS];  
    int count = 0;  

    // 加载数据到points数组中
    FILE *file = fopen(filename, "r");  

    while (fscanf(file, "%lf %lf", &points[count].x, &points[count].y) == 2) {  
        count++;  
        if (count >= MAX_POINTS) break;  
    }  
    fclose(file);  

// 计算参数  
     double sum_x = 0, sum_y = 0, sum_x2 = 0, sum_xy = 0;  
    
    for (int i = 0; i < count; i++) {  
        sum_x += points[i].x;  
        sum_y += points[i].y;  
        sum_x2 += points[i].x * points[i].x;  
        sum_xy += points[i].x * points[i].y;  
    }  

    double a = (count * sum_xy - sum_x * sum_y) / (count * sum_x2 - sum_x * sum_x);  
    double b = (sum_y - a * sum_x) / count;  

    printf("The fitted curve is y = %lfx + %lf\n", a, b);  
    // 计算残差平方和
    double sum1 = 0.0 , sum2 = 0.0 ,sum_loss = 0.0;
    for (int i = 0; i < count; i++) {
        double predicted_y = a * points[i].x + b;
        sum_loss += (points[i].y - predicted_y) * (points[i].y - predicted_y);
        sum1 += points[i].y / count;
        sum2 += predicted_y / count;
    }

    double t_loss = (sum1-sum2) * (sum1-sum2);
    double accuracy = (1 - sum_loss / t_loss) * 100;
    printf("Accuracy: %.1lf%% \n" , accuracy);
}
// 二次函数最小二乘法拟合  
void least_squares2(const char* filename) {  
    Point points[MAX_POINTS];  
    int count = 0;  

    // 加载数据  
    FILE *file = fopen(filename, "r");  


    while (fscanf(file, "%lf %lf", &points[count].x, &points[count].y) == 2) {  
        count++;  
        if (count >= MAX_POINTS) break;  
    }  
    fclose(file);  

    // 计算参数  
     double sum_x = 0, sum_y = 0, sum_x2 = 0, sum_x3 = 0, sum_x4 = 0, sum_xy = 0, sum_x2y = 0;  

    for (int i = 0; i < count; i++) {  
        sum_x += points[i].x;  
        sum_y += points[i].y;  
        sum_x2 += points[i].x * points[i].x;  
        sum_x3 += points[i].x * points[i].x * points[i].x;  
        sum_x4 += points[i].x * points[i].x * points[i].x * points[i].x;  
        sum_xy += points[i].x * points[i].y;  
        sum_x2y += points[i].x * points[i].x * points[i].y;  
    }  

     double D = count * (sum_x2 * sum_x4 - sum_x3 * sum_x3) - sum_x * (sum_x * sum_x4 - sum_x2 * sum_x3) + sum_x2 * (sum_x * sum_x3 - sum_x2 * count);  
     double a = (count * (sum_x2y * sum_x4 - sum_x3 * sum_y) - sum_x * (sum_xy * sum_x4 - sum_x2 * sum_y) + sum_x2 * (sum_xy * sum_x3 - sum_x2y * count)) / D;  
     double b = (count * (sum_xy * sum_x4 - sum_x3 * sum_x2y) - sum_x * (sum_xy * sum_x3 - sum_x2 * sum_y) + sum_x2 * (sum_x * sum_y - sum_x2y * count)) / D;  
     double c = (count * (sum_x2 * sum_x2y - sum_xy * sum_x3) - sum_y * (sum_x2 * sum_x4 - sum_x3 * sum_x2) + sum_x * (sum_xy * sum_x3 - sum_x2y * sum_x2)) / D;  

    printf("The fitted curve is y = %lfx^2 + %lfx + %lf\n", a, b, c);  
    double sum1 = 0.0 , sum2 = 0.0 ,sum_loss = 0.0;
    // 计算残差平方和
    for (int i = 0; i < count; i++) {
        double predicted_y = a * points[i].x * points[i].x + b * points[i].x + c;
        sum_loss += (points[i].y - predicted_y) * (points[i].y - predicted_y);
        sum1 += points[i].y / count;
        sum2 += predicted_y / count;
    }

    double t_loss = (sum1-sum2) * (sum1-sum2);
    double accuracy = (1 - sum_loss / t_loss) * 100;
    printf("Accuracy: %.1lf%% \n" , accuracy);
}

// 三次函数最小二乘法拟合  
void least_squares3(const char* filename) {  
    Point points[MAX_POINTS];  
    int count = 0;  

    // 加载数据  
    FILE *file = fopen(filename, "r");  


    while (fscanf(file, "%lf %lf", &points[count].x, &points[count].y) == 2) {  
        count++;  
        if (count >= MAX_POINTS) break;  
    }  
    fclose(file);  

    // 计算参数  
     double sum_x = 0, sum_y = 0, sum_x2 = 0, sum_x3 = 0, sum_x4 = 0, sum_x5 = 0, sum_x6 = 0,  
          sum_xy = 0, sum_x2y = 0, sum_x3y = 0;  

    for (int i = 0; i < count; i++) {  
        sum_x += points[i].x;  
        sum_y += points[i].y;  
        sum_x2 += points[i].x * points[i].x;  
        sum_x3 += points[i].x * points[i].x * points[i].x;  
        sum_x4 += points[i].x * points[i].x * points[i].x * points[i].x;  
        sum_x5 += points[i].x * points[i].x * points[i].x * points[i].x * points[i].x;  
        sum_x6 += points[i].x * points[i].x * points[i].x * points[i].x * points[i].x * points[i].x;  
        sum_xy += points[i].x * points[i].y;  
        sum_x2y += points[i].x * points[i].x * points[i].y;  
        sum_x3y += points[i].x * points[i].x * points[i].x * points[i].y;  
    }  

    // 构造方程组，解决以获得a, b, c, d  
     double D = count * (sum_x2 * (sum_x4 * sum_x6 - sum_x5 * sum_x5) - sum_x3 * (sum_x3 * sum_x6 - sum_x4 * sum_x5))  
              - sum_x * (sum_x * (sum_x4 * sum_x6 - sum_x5 * sum_x5) - sum_x2 * (sum_x3 * sum_x6 - sum_x4 * sum_x5))  
              + sum_x2 * (sum_x * (sum_x3 * sum_x5 - sum_x2 * sum_x4) - count * (sum_x3 * sum_x3 - sum_x2 * sum_x2));  

     double a = (count * (sum_x2y * (sum_x4 * sum_x6 - sum_x5 * sum_x5) - sum_xy * (sum_x3 * sum_x6 - sum_x4 * sum_x5))  
                - sum_x * (sum_xy * (sum_x4 * sum_x6 - sum_x5 * sum_x5) - sum_x2 * (sum_x2y * sum_x6 - sum_x3 * sum_x5))  
                + sum_x2 * (sum_x * (sum_xy * sum_x5 - sum_x2y * sum_x4) - count * (sum_xy * sum_x3 - sum_x2y * sum_x2))) / D;  

     double b = (count * (sum_xy * (sum_x4 * sum_x6 - sum_x5 * sum_x5) - sum_x2y * (sum_x3 * sum_x6 - sum_x4 * sum_x5))  
                - sum_x * (sum_xy * (sum_x3 * sum_x6 - sum_x4 * sum_x5) - count * (sum_x2y * sum_x6 - sum_x3 * sum_x5))  
                + sum_x2 * (count * (sum_x2y * sum_x5 - sum_xy * sum_x4) - sum_xy * (sum_xy * sum_x3 - sum_x2y * sum_x2))) / D;  

     double c = (count * (sum_x2 * (sum_x2y * sum_x4 - sum_xy * sum_x5)) - sum_y * (sum_x2 * sum_x4 - sum_x3 * sum_x3)) / D;  
     double d = (sum_y - a * sum_x2 - b * sum_x - c * count) / count;  // 根据y的总和得到d  

    printf("The fitted curve is y = %lfx^3 + %lfx^2 + %lfx + %lf\n", a, b, c, d);  
    double sum1 = 0.0 , sum2 = 0.0 ,sum_loss = 0.0;
    // 计算残差平方和
    for (int i = 0; i < count; i++) {
        double predicted_y = a * points[i].x * points[i].x * points[i].x + b * points[i].x * points[i].x + c * points[i].x + d;
        sum_loss += (points[i].y - predicted_y) * (points[i].y - predicted_y);
        sum1 += points[i].y / count;
        sum2 += predicted_y / count;
    }

    double t_loss = (sum1-sum2) * (sum1-sum2);
    double accuracy = (1 - sum_loss / t_loss) * 100;
    printf("Accuracy: %.1lf%% \n" , accuracy);
}
void gradient_descent(const char* filename) 
{
    Point points[MAX_POINTS];  
    int count = 0;  

    // 加载数据到points数组中
    FILE *file = fopen(filename, "r");  

    while (fscanf(file, "%lf %lf", &points[count].x, &points[count].y) == 2) {  
        count++;  
        if (count >= MAX_POINTS) break;  
    }  
    fclose(file);  

    double a = 0;//初始化a 
    double b = 0;//初始化b,初始化为合理范围内的实数即可，取得的结果稍有不同，基本可以认为是一致的 
	//不断调整迭代a、b大小，直至函数走到最低点，解出y=ax+b 
    double lr = 0.001;//初始化学习率，即梯度下降的步长 
    
    //梯度下降部分 
    for (int epoch = 0; epoch < EPOCHS; epoch++) {
        double dJ_da = 0;//对a求偏导，初始化 
        double dJ_db = 0;//对b求偏导，初始化 
        for (int i = 0; i < count; i++) {
            dJ_da -= 2 * points[i].x * (points[i].y-(a * points[i].x+b));//对a求偏导(还没除以N)
            dJ_db -= 2 * (points[i].y-(a * points[i].x+b));//对b求偏导(还没除以N)
            //printf("%f %f\n",dJ_da,dJ_db);注释掉的这行检验拟合状态 
			//如果学习率和最大迭代次数设置恰当，则最终偏导数应该是处于时正时负的摇摆状态
			//即认为最后阶段在最低点左右徘徊，梯度下降经常用下坡来打比方
			//迭代的最终目的是找到最低点，即在损失函数取最小时得到a、b、c的值 
        }
        dJ_da /= count;
        dJ_db /= count;//再循环外做最后一次处理，得到偏导，此处不能用int 
        a -= lr * dJ_da;//步长乘以偏导为在a轴移动的距离 
        b -= lr * dJ_db;//在b轴移动的距离 
    }//循环，得到a，b
	
	printf("The fitted curve is y = %lfx + %lf\n", a, b); 
    // 计算残差平方和
    double sum1 = 0.0 , sum2 = 0.0 ,sum_loss = 0.0;
    for (int i = 0; i < count; i++) {
        double predicted_y = a * points[i].x + b;
        sum_loss += (points[i].y - predicted_y) * (points[i].y - predicted_y);
        sum1 += points[i].y / count;
        sum2 += predicted_y / count;
    }

    double t_loss = (sum1-sum2) * (sum1-sum2);
    double accuracy = (1 - sum_loss / t_loss) * 100;
    printf("Accuracy: %.1lf%% \n" , accuracy);

}

void gradient_descent2(const char* filename) 
{
    Point points[MAX_POINTS];  
    int count = 0;  

    // 加载数据到points数组中
    FILE *file = fopen(filename, "r");  

    while (fscanf(file, "%lf %lf", &points[count].x, &points[count].y) == 2) {  
        count++;  
        if (count >= MAX_POINTS) break;  
    }  
    fclose(file);  
    
    // 不断调整迭代a、b、c的大小，直至函数走到最低点，解出y=ax^2+bx+c
    double a = 5; // 初始化a
    double b = 0; // 初始化b
    double c = 0; // 初始化c
    double lr = 0.001; // 初始化学习率，即梯度下降的步长
    //__mingw_printf("%.12Lf",x[0]); 在重复输出0后尝试寻找原因，发现是编译器的问题 

    // 梯度下降部分
    int max_iterations = 1000; // 最大迭代次数，达到最大迭代次数时停止迭代
    for (int iter = 0; iter < max_iterations; iter++) {
        double dJ_da = 0; // 对a求偏导，初始化
        double dJ_db = 0; // 对b求偏导，初始化
        double dJ_dc = 0; // 对c求偏导，初始化
        
        for (int i = 0; i < count; i++) {
            dJ_da -= 2.0 * points[i].x * points[i].x * (points[i].y - (a * points[i].x * points[i].x + b * points[i].x + c)); // 对a求偏导(还没除以N)
            dJ_db -= 2.0 * points[i].x * (points[i].y - (a * points[i].x * points[i].x + b * points[i].x + c)); // 对b求偏导(还没除以N)
            dJ_dc -= 2.0 * (points[i].y - (a * points[i].x * points[i].x + b * points[i].x + c)); // 对c求偏导
            //这里利用了分别求导可以累加的特点简化原方程的运算 
           // __mingw_printf("%Lf %Lf %Lf\n",dJ_da,dJ_db,dJ_dc);
        }
        dJ_da /= count;
        dJ_db /= count; // 在循环外做最后一次处理，得到偏导，此处不能用int
        dJ_dc /= count;
        a -= lr * dJ_da; // 步长乘以偏导为在a轴移动的距离
        b -= lr * dJ_db; // 在b轴移动的距离
        c -= lr * dJ_dc; // 在c轴移动的距离
    } // 循环，得到a，b，c

    printf("The fitted curve is y = %lf x^2 + %lf x + %lf\n", a, b, c);
    double sum1 = 0.0 , sum2 = 0.0 ,sum_loss = 0.0;
    // 计算残差平方和
    for (int i = 0; i < count; i++) {
        double predicted_y = a * points[i].x * points[i].x + b * points[i].x + c;
        sum_loss += (points[i].y - predicted_y) * (points[i].y - predicted_y);
        sum1 += points[i].y / count;
        sum2 += predicted_y / count;
    }

    double t_loss = (sum1-sum2) * (sum1-sum2);
    double accuracy = (1 - sum_loss / t_loss) * 100;
    printf("Accuracy: %.1lf%% \n" , accuracy);

}


// 三次函数梯度下降法拟合  
void gradient_descent3(const char* filename) {  
    Point points[MAX_POINTS];  
    int count = 0;  

    // 加载数据  
    FILE *file = fopen(filename, "r");  


    while (fscanf(file, "%lf %lf", &points[count].x, &points[count].y) == 2) {  
        count++;  
        if (count >= MAX_POINTS) break;  
    }  
    fclose(file);  

    // 初始化参数  
     double a = 0, b = 0, c = 0, d = 0; // 参数 a, b, c, d  
     double learning_rate = 0.001; // 学习率  
    int iterations = 10000; // 迭代次数  

    // 梯度下降迭代计算参数  
    for (int iter = 0; iter < iterations; iter++) {  
         double gradient_a = 0;  
         double gradient_b = 0;  
         double gradient_c = 0;  
         double gradient_d = 0;  

        for (int i = 0; i < count; i++) {  
             double x = points[i].x;  
             double y = points[i].y;  
             double prediction = (a * x * x * x) + (b * x * x) + (c * x) + d;  
             double error = prediction - y;  

            gradient_a += (3.0 / count) * error * x * x; // 对a的梯度  
            gradient_b += (2.0 / count) * error * x; // 对b的梯度  
            gradient_c += (1.0 / count) * error; // 对c的梯度  
            gradient_d += (1.0 / count) * error; // 对d的梯度  
        }  

        a -= learning_rate * gradient_a; // 更新a  
        b -= learning_rate * gradient_b; // 更新b  
        c -= learning_rate * gradient_c; // 更新c  
        d -= learning_rate * gradient_d; // 更新d  
    }  

    printf("The fitted curve is y = %Lfx^3 + %lfx^2 + %lfx + %lf\n", a, b, c, d);  
    double sum1 = 0.0 , sum2 = 0.0 ,sum_loss = 0.0;
    // 计算残差平方和
    for (int i = 0; i < count; i++) {
        double predicted_y = a * points[i].x * points[i].x * points[i].x + b * points[i].x * points[i].x + c * points[i].x + d;
        sum_loss += (points[i].y - predicted_y) * (points[i].y - predicted_y);
        sum1 += points[i].y / count;
        sum2 += predicted_y / count;
    }

    double t_loss = (sum1-sum2) * (sum1-sum2);
    double accuracy = (1 - sum_loss / t_loss) * 100;
    printf("Accuracy: %.1lf%% \n" , accuracy);
}  




// 生成数据集部分  
void gendat(int t,   double a,     double b,   double c,     double d, const char* filename)   
{  
    FILE *file = fopen(filename, "w");  
    // if (!file) {  
    //     perror("Error opening file");  
    //     return;  
    // }  
    srand(time(NULL)); // 初始化随机数种子  

    for (int i = 0; i < MAX_POINTS; i++) {  
         double x =      (double)rand() / RAND_MAX * 200 - 100; // 生成-100到100之间的随机数  
         double y = 0;  

        switch (t) {  
            case 1: // 一次函数  
                y = a * x + b;  
                break;  
            case 2: // 二次函数  
                y = a * x * x + b * x + c;  
                break;  
            case 3: // 三次函数  
                y = a * x * x * x + b * x * x + c * x + d;  
                break;  
        }  

        // 添加噪声  
        y += (rand() / (RAND_MAX / 10.0)) - 5;  

        fprintf(file, "%lf %lf\n", x, y); // 写入文件  
    }  

    fclose(file);  
}  




int main() 
{  
    int t;  
    double a, b, c, d;  
    char method[10];  
    
    printf("function type : (1, 2, 3): ");  //指定函数类型为最高项t次
    scanf("%d", &t);  
    printf("parameters : a, b, c, d: ");  //参数 没有输入0
    scanf("%lf %lf %lf %lf", &a, &b, &c, &d);  
    printf("method : (ls or gd): ");  //拟合方法 ls:最小二乘法 gd:梯度下降法 
    scanf("%s", method);  

    // 生成数据集  
    gendat(t, a, b, c, d, "some.dat");  

    //选择拟合方式  
    if (strcmp(method, "ls") == 0) {  
        if (t == 1) least_squares("some.dat"); 
        else if (t == 2)  least_squares2("some.dat"); 
        else least_squares3("some.dat");
    } 
    else if (strcmp(method, "gd") == 0) {  
        if (t == 1) gradient_descent("some.dat"); 
        else if (t == 2)  gradient_descent2("some.dat"); 
        else gradient_descent3("some.dat");
    } 

    return 0;  
}